The Fundamental Matrix of the System of Linear Elastodynamics in Hexagonal Media. Solution to the Problem of Conical Refraction

An explicit integral representation by single definite integrals of the fundamental matrix (Green’s tensor) of the time-dependent system of hexagonal elastic media is derived. Thereby the problem of internal conical refraction in such media is solved. 0. Introduction and notations Hexagonal or transversely isotropic media are characterized by the property of rotational symmetry with respect to an axis. In this paper, we extend R. G. Payton’s seminal work in that area by providing qualitative and quantitative information on the fundamental matrix E (= “free space Green’s tensor”) of the elastodynamic system P (∂) = P (∂t,∇) = I3∂ t +A(∇) for hexagonal crystals. On the one hand, we determine the singular support of E (Prop. 4), thereby solving the open problem of conical refraction; on the other hand, we represent E by simple definite integrals (Prop. 5). Up to now, representations of E were given by loop integrals over curves defined as intersections of the slowness surface detP (1, ξ) = 0 with the planes Π(t,x) = {ξ ∈ R; t + xξ = 0}, see [5, (5.4), p. 45], [31, (6.7), (6.8), p. 3318], [38, (A), p. 327], or formula (8) in Prop. 3 below. Though these loop integrals running over curves of genus 3 (in general) yield some geometric insight into the structure of E, they are harder to evaluate numerically than the simple integrals in Prop. 5. Also, the specialization to particular cases is not as straightforward, cf. e.g. [38, Section 4, p. 336]. 1 2 N. ORTNER AND P. WAGNER In the recent article [7], R. Burridge and J. Qian take the loop integral approach for the system of crystal optics. Let us mention, incidentally, that the “static term” (see [7, (5.11), p. 78]) was calculated for the magnetic field in [38, p. 335], cf. also [49, p. 415]. We also remark that, for the electric field, the static term E1 can be represented more simply than in [7], namely E1 = − Y (t)t 4πdetσ ∇∇ T 1 √ xTσ−1x . Our representation of E starts from a variant of the Herglotz–Petrovsky–Leray formula for the fundamental matrix of a system with a strictly hyperbolic determinant, which formula we call the Herglotz–G̊arding formula (see § 2, Prop. 1). In our analysis of the qualitative properties of E in § 3, we first give necessary and sufficient conditions for the hyperbolicity of the operator P (∂) in terms of the elastic constants (see Prop. 2), and then adjust the formula of Prop. 1 to the system of hexagonal elastodynamics, which is not strictly hyperbolic (Prop. 3). Next we determine in Prop. 4 the singular support of E, i.e. the set where E is not smooth. Its intersection with t = 1, the so-called wave-front surface is, essentially, the dual of the slowness surface. However, “conical” (i.e. singular) points on the slowness surface can result in (internal) conical refraction. Mathematically, this amounts to additional plane or conical lids contained in the wave-front surface. (These lids are given by the supports of the corresponding localized operators.) For scalar operators in three space dimensions, these lids are always present due to the foundational results of M. F. Atiyah, R. Bott, and L. G̊arding (cf. e.g. [2, Th. 7.7, p. 175]), whereas for systems, these lids can disappear. (This phenomenon is connected with the form of the adjoint matrix of P ; it is most easily understood when considering a diagonal 2× 2 system P = ( P1 0 0 P2 ) : Points ξ0 which belong to the intersection of the slowness surfaces P1(1, ξ) = 0, P2(1, ξ) = 0 and are regular on both of them are singular on the slowness surface of P (∂), i.e. of the determinant operator detP (∂) = P1(∂)P2(∂), without leading to conical refraction for the system P (∂). In fact, P (∂) has the fundamental matrix ( E1 0 0 E2 ) if Ei is the fundamental solution of Pi(∂), i = 1, 2.) For a general account on the history and mathematics of conical refraction, cf. [17] and [35]. We shall show that conical refraction does not occur in the propagation of elastic waves in general hexagonal media (see Prop. 4 and the Remark 1 thereupon). This qualitative analysis is related to the historical development as follows: In 1954 and 1957, J. de Klerk and M. J. P. Musgrave predicted conical refraction for transverse elastic waves in cubic and in tetragonal crystals as a result of their investigation of the conical points on the slowness surface. A proof of the existence of conical refraction was given for cubic crystals by R. Burridge [5] in 1967, and for tetragonal crystals by P. A. Barry and M. J. P. Musgrave in 1978 [3]. In 1983, R. G. Payton [41, pp. 66, 67] conjectured that conical refraction did not occur in transversely isotropic elastic solids. However, in 1992 [42], he proved the existence of conical refraction in special hexagonal media fulfilling c11 = c33 = c44. The peculiarity of the cases c11 = c44, c33 = c44, and c13 + c44 = 0, respectively, was also recognized and investigated by P. Chadwick and A. L. Shuvalov [11] in 1997. As mentioned above, we shall show that no wave-front lids exist in general transFUNDAMENTAL MATRICES OF HEXAGONAL MEDIA 3 versely isotropic solids, and we shall characterize precisely the limit cases in which such lids (i.e. conical refraction) appear. The definite integrals we obtain in Prop. 5 for the fundamental matrix are, in general, complete Abelian integrals of genus 3. The calculation of the genus of the complex algebraic curve connected with the integral representation (cf. Rem. 1 to Prop. 5) clarifies the nature of these integrals. They reduce to integrals of genus 0 (i.e. algebraic functions) iff c13 + c44 = 0 or (c13 + c44) 2 = (c11 − c44)(c33 − c44), for which cases we deduce in § 5 the known results of R. G. Payton ([41], cf. also [6], [38]). Furthermore, a reduction of the general representation formula to elliptic integrals (i.e. genus 1) occurs if c11 = c44, a particular case of which was studied in [42]. In § 6, the values of the fundamental matrix in the x1, x2 -plane are expressed by elliptic integrals. Different representation formulas for this case were already deduced in 1967 by N. Cameron and G. Eason [8]. Finally, we calculate the values of E on the x3 -axis (see Prop. 6) and obtain the results of R. G. Payton [40], [41]. We make constant use of L. Schwartz’ theory of distributions ([46]), since we agree with L. G̊arding that, before the invention of distribution theory, “the analysis of singularities of solutions to partial differential equations was a painful matter, especially for non-elliptic ones. A real understanding first became possible with the framework of Schwartz’s distributions ...” ([15, p. 32]). Finally, let us establish some notations. As usual, Euclidean space is written as R and Sn−1 denotes the unit sphere in R, i.e., Sn−1 := {x ∈ R : |x| = 1}. In general, the variables in R are written as (t, x) = (t, x1, . . . , xn−1) or (τ, ξ) = (τ, ξ1, . . . , ξn−1). The Heaviside function is denoted by Y, the l × l unit matrix is written as Il, and A ad is the adjoint matrix of the l × l matrix A, i.e., A · A = (detA) · Il. We consider differential operators with constant coefficients only and use as differentiation symbols ∂t := ∂ ∂t , ∇ := (∂1, . . . , ∂n−1) := ( ∂ ∂x1 , . . . , ∂ ∂xn−1 ) , ∂ := (∂t,∇), ∆l := ∂ 1 + · · ·+ ∂ l . We employ the standard notations for the distribution spaces D′, S ′, the duals of the spaces D, S of “test functions” and of “rapidly decreasing functions”, respectively, cf. [13], [19], [26], [28], [45], [46], and we use angle brackets, e.g. ⟨φ, T ⟩, to indicate the evaluation of the distribution T on the test function φ. By suppT, sing suppT, sing suppA T, we denote the support, the singular support, and the analytic singular support of T, respectively, i.e. the complements of the sets, where T vanishes, is infinitely differentiable, or is real-analytic, respectively, cf. [26, 2.2 and 8.4], [13, 1.4.1, 8.6.1]. We use the Fourier transform F in the form (Fφ)(x) := ∫ e−ixξφ(ξ) dξ, φ ∈ S, this being extended to S ′ by continuity. Herein and also elsewhere, the Euclidean inner product (x, ξ) 7→ xξ is simply expressed by juxtaposition. 4 N. ORTNER AND P. WAGNER In several places, we use the abbreviations ρ := ξ 1 + ξ 2 2 , x ′ := (x1, x2), H+ := {(t, x) ∈ R; t ≥ 0} and Π(t,x) := {ξ ∈ R; t+ xξ = 0}. For the convenience of the reader, the constants a1, . . . , a7 and the wave operators W1, . . . ,W5, which appear scattered in the text, are collected in an appendix. 1. The equations of linear elastodynamics in transversely isotropic media Let us recall the equations governing the displacement in transversely isotropic media, cf. [41, (1.1.6), (1.3.2), pp. 1, 3], [35, (4.5.1b), p. 57, p. 94]. In anisotropic media, the laws of Newton and of Hooke yield (1) ( ρI3∂ 2 t +A(∇) ) u = ρf, where ρ denotes the (constant) density of mass, u and f represent the vectors of displacement and of exterior forces (per unit mass), respectively, and A(∇) is a matrix of second-order linear differential operators determined by the tensor (cpqrs) of the elastic constants. The dimension of the linear space of tensors (cpqrs) fulfilling the appropriate symmetry relations equals 21 in general, but reduces to 5 in the case of transversely isotropic media. Following [4, (6.7), (6.10), pp. 572, 573] and [38, 1.2.3, p. 319], we abbreviate (2) a1 = c1111 = c11, a2 = c3333 = c33, a3 = c1133 + c2323 = c13 + c44, a4 = 1 2 (c1111 − c1122) = 1 2 (c11 − c12), a5 = c2323 = c44, (cij being the “contracted index notation”, cf. [41, (1.3.1), p. 3]), and then obtain (3) A(∇) = −  a1∂ 1 + a4∂ 2 + a5∂ 3 (a1 − a4)∂1∂2 a3∂1∂3 (a1 − a4)∂1∂2 a4∂ 1 + a1∂ 2 + a5∂ 3 a3∂2∂3 a3∂1∂3 a3∂2∂3 a5(∂ 2 1 + ∂ 2 2) + a2∂ 2 3  . Without restriction of generality, we set the mass density ρ equal to 1. The solution u of (1) in unbounded space is given as a convolution integral of the fundamental matrix E with f. Here and in the following, E denotes the uniquely determined matrix of distributions with support in the half-space H+ := {(t, x) ∈ R; t ≥ 0} fulfilling the matrix equation ( I3∂ 2 t +A(∇) ) E = I3δ, where δ denotes the Dirac measure at the origin. Hence the i -th column of E corresponds to the exterior force vector eiδ with I3 = (e1|e2|e3). (Note that there is no generally accepted term for the concept of fundamental matrix. In physics, it is often called “Green’s tensor” (cf. [41]) or “free space Green’s tensor”, in the mathematical literature “noyau élémentaire à droite” (cf. [46, p. 140]), “fundamental solution to the right” (cf. [25, 3.8, p. 94]), “Green’s matrix” (cf. [20, III.1, p. 106]), “fundamental solution” (cf. [16, p. 215]).) FUNDAMENTAL MATRICES OF HEXAGONAL MEDIA 5 Similarly as in [12], [9], [35], [41], we shall illustrate our calculations by two specific transversely isotropic materials, namely cobalt and titanium boride (TiB2), which are characterized by the following elastic constants (in gigapascal): cobalt: a1 = 307, a2 = 358, a3 = 178.5, a4 = 71, a5 = 75.5 TiB2 : a1 = 690, a2 = 440, a3 = 570, a4 = 140, a5 = 250 see [41, Table 1, p. 3], [9, Table 4, p. 37]. 2. The Herglotz–G̊arding formula for the fundamental matrix of hyperbolic systems For the explicit representation of the fundamental matrix of a homogeneous system P (∂) with strictly hyperbolic determinant, we shall deduce a modification of the Herglotz–Petrovsky–Leray formula (see [1, pp. 173–177], [27, (12.6.10), p. 129]), which modification we call the Herglotz–G̊arding formula, since it is derived from the formula in [17, Thm. 2, p. 375] for the fundamental solution of a homogeneous, strictly hyperbolic scalar operator. P (∂) denotes an l × l matrix of constant coefficient differential operators acting on R and homogeneous of degree m and we always set Q := detP. (In § 1, P (∂) corresponds to I3∂ 2 t + A(∇).) P (∂) is assumed to be hyperbolic with respect to t, i.e. (i) Q(1, 0) ̸= 0, and (ii) the polynomial τ 7→ Q(τ, ξ) has only real roots for each ξ ∈ Rn−1, cf. [1, p. 129], [29, pp. 89, 90], [27, Thm. 12.4.3, p. 113]. If these roots are pairwise different for each ξ ∈ Rn−1 \ {0}, then Q is called strictly hyperbolic, cf. [1, Def. 3.8, p. 129], [27, Def. 12.4.11, p. 118]. (Note that the strict hyperbolicity of P, i.e. the hyperbolicity of P + P1 for all matrices P1 of differential operators of lower order, is implied by, but is not equivalent to, the strict hyperbolicity of Q, cf. [30, p. 787], [16, p. 221].) Proposition 1. Let P (τ, ξ) = P (τ, ξ1, . . . , ξn−1) be a real l× l matrix of polynomials which are homogeneous of degree m and suppose that Q(∂) = detP (∂) is strictly hyperbolic with respect to t and that Q(τ, ξ) does not contain τ as a factor. Define the measure T ∈ D′(Rn \ {0})l×l by T := P (τ, ξ) δ ( Q(τ, ξ) ) sign ( (∂τQ)(τ, ξ) ) . Furthermore, set s+ := Y (s)s λ ∈ Lloc(Rs) for Reλ > −1 and let s n−m−1 + be the finite part evaluated at n−m− 1 of the meromorphic extension to the whole complex plane of the holomorphic function {λ ∈ C; Reλ > −1} −→ S ′(R) : λ 7−→ s+, 6 N. ORTNER AND P. WAGNER cf. [19, pp. 48, 49]. Then the uniquely determined fundamental matrix E of P (∂) with support in H+ fulfills (4) E(t, x) = −2(2π)1−nY (t) ∫ Rn−1 T (1, ξ)Re [ im+1Fsn−m−1 + ] (t+ xξ) dξ+ Y (t)S(t, x), where S = 0 if n is even or m < n, and is otherwise an l×l matrix of homogeneous polynomials of degree m− n. Remarks. 1) Due to the homogeneity of T, the restriction of T to the hyperplane τ = 1 is well-defined; furthermore, the integral ∫ Rn−1 · · · dξ in (4) has to be understood in the distributional sense, i.e. for φ ∈ D(Rt,x) with suppφ ⊂ H+, we have ⟨φ,E⟩ = −2(2π)1−n ∫ Rn−1 T (1, ξ)ψ(ξ) dξ + ⟨φ, S⟩ ∈ Cl×l where ψ(ξ) := ⟨φ(t, x),Re [ im+1Fsn−m−1 + ] (t+ xξ)⟩ ∈ C∞(Rn−1 ξ ) and T (1, ξ)ψ(ξ) is an l × l matrix of integrable measures. Note also that the multiplication with Y (t) in formula (4) is well-defined, since the support of the following distribution intersects the hyperplane t = 0 in the origin x = 0 only, and since a homogeneous distribution of degree m − n can uniquely be continued from R \ {0} to R (for m ≥ 1). 2) Formula (4), which we call the Herglotz–G̊arding formula, generalizes earlier formulas such as [23, (4)–(13), pp. 609, 610], [24, (7.58), p. 192], [1, pp. 176, 177], [17, Thm. 2, p. 375], [38, Thm., p. 324], [51, Prop. 1, p. 309]. Proof. a) By the hyperbolicity of P, the matrices P (τ ± iε, ξ) are invertible for ε > 0 and (τ, ξ) ∈ R, and the entries of the inverse matrices P (τ ± iε, ξ)−1 grow at most polynomially when ε ↘ 0. Hence the two limits limε↘0 P (τ ± iε, ξ)−1 exist in D′(Rn)l×l (cf. [1, p. 121]) and yield homogeneous distributions of degree −m. Since F−1 ( lim ε↘0 P (τ ± iε, ξ)−1 ) are the two fundamental matrices of P (−i∂t,−i∇) with support in ∓H+ respectively (cf. [27, (12.5.3), p. 120]), we obtain E = Y (t)i1−m2πF−1T with T = 1 2πi lim ε↘0 ( P (τ − iε, ξ)−1 − P (τ + iε, ξ)−1 ) ∈ S ′(Rn)l×l. Next, due to the strict hyperbolicity of Q(∂) = detP (∂), Sokhotsky’s formula lim ε↘0 1 x± iε = vp 1 x ∓ iπδ in D′(R1) FUNDAMENTAL MATRICES OF HEXAGONAL MEDIA 7 implies (outside the origin) T = 1 2πi P (τ, ξ) lim ε↘0 ( Q(τ − iε, ξ)−1 −Q(τ + iε, ξ)−1 ) = P (τ, ξ) δ ( Q(τ, ξ) ) sign(∂τQ)(τ, ξ), cf. [38, p. 322]. Note that the last expression is defined in D′(Rn \ {0}) in the usual way: ⟨φ, δ ( Q(τ, ξ) ) sign(∂τQ)(τ, ξ)⟩ =